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Abstract 

Droplet interaction with a high temperature gaseous crossflow is important because of its wide 
application in systems involving two phase mixing such as in combustion requiring quick mixing 
of fuel and air. The focus of this work is to investigate dispersion of a two-dimensional 
evaporating spray into a crossflow. 

An interactive Microsoft® Excel program for tracking a single droplet in crossflow that has 
previously been developed was modified to include droplet evaporation computation. In addition 
to the high velocity airflow, the injected droplets are also subjected to increased combustor 
temperature and pressure that affect their motion in the flow field. Six ordinary differential 
equations (namely the time rate of change of x, z, Ud, wj, D, and T s ) are then solved by 4 th -order 
Runge-Kutta method using Microsoft® Excel software. 

Visual Basic programming and Excel macrocode are used to calculate the data and plot the 
droplet’s motion in the flow field. This program computes and plots the data sequentially without 
forcing the user to open other types of plotting programs. A user’s manual on how to use the 
program is included. 
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Symbol List 


Ad Projected area of the droplet 

B m Mass transfer number 

Bt Heat transfer number 

C exp Expansion coefficient 

C D Drag coefficient 

c p Specific heat at constant pressure 

c p ,a Specific heat of air 

c Pj f Specific heat of liquid fuel 

c p ,fv Specific heat of fuel vapor 

c p ,g Specific heat of fuel-air mixture 

D Droplet diameter 

D F a Diffusion coefficient of fuel in air 

Do Initial diameter of droplet 

F Force 

g Gravitational acceleration constant 

kA Thermal conductivity of air 

k Fv Thermal conductivity of fuel vapor 

k g Thermal conductivity of fuel-air mixture 

F Fatent heat of vaporization 

Fxbn Fatent heat at nonnal boiling point 

M a Molecular weight of air 

Mf Molecular weight of fuel 

nidrop Mass of a droplet 

m Evaporation rate 

m " Evaporation rate per unit surface area 

n Exponential constant 

Nu Nusselt number 

P Ambient pressure 

Pf Partial pressure of fuel 

P va P Vapor pressure 

PrA Prandtl number of air 

Qact Total heat transfer to droplet 

Q ss Steady-state heat transfer to droplet 

R u Universal gas constant 

Red Reynolds number of droplet 

r Spherical radius of the droplet 

r 32 0.5 *(Sauter mean radius) 

Tbn Nonnal boiling point 
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T cr it Critical temperature 

Tdrop Droplet temperature 

T r Reference temperature 

T s Droplet surface temperature 

Too Ambient temperature 

t e Evaporation time 

U R Relative velocity between the droplet and the gas stream 
ua Velocity of the cross stream (air) in the x-direction 
Ud Velocity of the droplet in the x-direction 
Vd Droplet Volume 

wa Velocity of the cross stream (air) in the z-direction 
Wd Velocity of the droplet in the z-direction 

x Horizontal direction (+ equals to the right) 

Ya Mass fraction of air 

Yp Mass fraction of fuel 

Z Vertical direction (+ equals up) 

a g Thermal diffusivity of gas 

A, Evaporation constant 

p A Viscosity constant of air 

Pa Density of the cross stream (air) 

pd Density of the droplet 

p F28S6K Density of fuel droplet at 288. 6K 
At time step-size 

Subscripts 

a Air 

d Droplet 

f Fuel 

g Mixture of gaseous phase 

r Reference condition 

s Droplet surface 

st Steady state 


NASA/TM— 2004-212910 


3 



1. Introduction 


A liquid spray injected into a gaseous crossflow with high temperature is important because of 
its wide application in systems involving two phase mixing. It is important to be able to compute 

this flow to optimize the mixing strategy. 

An existing Excel program (ref. 1) has previously been developed for tracking a single droplet in 
crossflow computation. This work is focused on producing a quick computational method for 
determining spray penetration with evaporation. With this spreadsheet, one can investigate the 
dispersion of an air-blast atomized spray jet into a high temperature crossflow. During the 
transverse injection of a spray into high velocity airflow, the droplets (carried along in the 
gaseous stream of co-flowing air) are not only subjected to forces due to the crossflow motion, 
but also to increases in the combustor temperature and pressure (see fig. 1). 



Figure 1.1 (from NASA/CR— 2000-210467) (ref. 2) 


2. Governing Equations 

2.1 Droplet Trajectories and Velocities 

The trajectories of the droplets can be tracked by applying a Lagrangian-based analysis to the 
droplets. The momentum equations for a droplet can be obtained by equating the droplet motion 
to: (1) The viscosity and pressure-related drag forces, and (2) The pressure gradient and viscous 
forces related to the fluid surrounding the droplet, and (3) The inertia of the virtual mass, induced 
when the particle acceleration affects the fluid mass acceleration. 

Droplet trajectory and velocity with respect to time can be calculated. Based on these principles 
along with the following assumptions: 

(1) The droplets are spherical, 

(2) No droplet breakup occurs, 

(3) Vaporization is not considered yet and will be derived separately in section 2.2, 

(4) Lift, virtual mass, and Basset forces which takes into account the acceleration history of 
the droplet, are neglected, 

(5) Chemical reaction is not included, 
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These assumptions reduce the droplet momentum equation to include only the effects of the drag 
and body forces. The general momentum equations for a single droplet injected along the 
positive x-direction, transversely into a downward-flowing air stream in the positive z-direction, 
as shown in figure 2, is described by 


F(l ~ F drag 


+ F, 


body 


( 1 ) 


where the net force F d that drives the droplet motion is balanced by the drag force opposing its 
motion, and the field forces acting on the droplet. The aerodynamic drag force is given by 


F drag 2 PgU R U h 


A d C D 


( 2 ) 


where p A is the air density, and Ad and Co, the projected area and the drag coefficient of the 
droplet, respectively. The relative velocity between the droplet and the crossflow has a 
magnitude of Ur (see fig. 2). The subscript “c/” refers to the droplet and “g” the crossflow air. 



The body force, resulting from an equivalent volume of air that buoys the droplet, includes the 
gravitational and buoyancy forces. It is given by 

^body = iPd ~ P/i ) Ki S (3) 

which says that the body force is equal to the product of relative droplet and air density (pj - p A ), 
the droplet volume V A , and the gravitational acceleration g. 
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Substituting equation (2) and equation (3) to equation (1) yields: 


Pd v d /. = U A ) U r 

at 2 

A d C o 

(4) 

Pd V d / = r.pj^d W a) U R A d C D + 
at 2 

{p A - pd) y d g 

(5) 

s 

II 


(6) 

dz 

— = w d 
dt 


(7) 


The drag coefficient of the droplet depends on the droplet Reynolds number and is given by 


24 , 1 2/3 

1 + -Re rf 

Re.L 6 _ 

Re rf < 1000 

(8) 

0.424 

Re rf > 1000 



where Rej is the droplet Reynolds number and is defined as follows 


Re^ =■ 


2 P. 


U t 


Pa 


(9) 


in which r d is the droplet radius and ju g is the gas (air) viscosity. 


2.2 Droplet Evaporation 

To include the effect of evaporation rate on spray penetration, apply a control volume at droplet 
surface that will change with droplet radius during evaporation process. For simplicity, consider 
steady state analysis first. 


2.2.1. Steady-State Analysis 

A fuel droplet rarely reaches a steady state (ref. 3) during its lifetime. This is because most 
commercial fuels are multi-component, where different fuel compounds posses its own 
properties, for example kerosene and gasoline. To simplify analysis, ‘steady state’ tenn here 
refers to ‘quasi-steady’, which allows droplet lifetime and evaporation rate to be estimated to an 
acceptable level of accuracy. 
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To simplify the analysis, in addition to the assumptions listed in trajectory analysis as shown in 
section 2.1: 

1 . there is no radiation, 

2. there is no internal circulation and internal convective heating within droplet, 

3. consider only single-component fuel (with well-defined boiling point), and 

4. it is quasi-steady flow. 

Consider a fuel droplet at low fuel injection temperature that is suddenly exposed to a gaseous 
crossflow at high temperature. Initially, almost all heat supplied to the droplet serves to raise the 
droplet temperature. As the droplet temperature rises, fuel vapor will fonn at the droplet surface 
and has two main effects: 

1 . a large portion of heat transferred to droplet is used to vaporize the droplet, and 

2. the outward flow of fuel vapor impedes the rate of heat transfer to droplet 

Eventually, a stage is reached where all heat transferred to droplet is used as the heat of 
vaporization and the droplet temperature will stabilize at a steady-state temperature. 


Mass Transfer Number 


Assume that the thermal diffusion is negligible. Therefore, the concentration gradient is the only 
driving force considered for species diffusion in the direction of the diffusion path. Then, the 
following expression for an evaporating fuel droplet of radius r is described by: 


dY F 

dr 


■^(m F Y A ) 

n p 

1J FA r 


( 10 ) 


where D FA 
m F 

P 

R u 

T 


Y,(r) 

Y a 


diffusion coefficient of fuel in air 

mass rate of diffusion per unit area (mass flux) 

ambient air pressure 

universal gas constant 

ambient air temperature 

fuel mass fraction 

air mass fraction at range r s <r< at any time 


Y a = 1-Y f 


(10a) 


From the continuity equation applied on the control surface surrounding a droplet, one obtains 


m, = m F 


f r \ 2 


\r ) 


( 11 ) 
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Figure 2.2 - Sketch of the control surface 
surrounding a droplet 

where m" s - mass flux at droplet surface 

r s - radius of droplet 
r - radius of control surface at time t 


Substituting equation (10a) and (11) into (10) yields: 


dY, 

dr 



where Y f s - fuel mass fraction at droplet surface 


Assume the ideal gas relation ( p = ), separating variables, integrating, and 

R u T 

equation (12) yields: 


m 


rr 

F,s 



•ln(l-F Fl ) 


where p - fuel density 

Multiplying by droplet surface area (A s = 4m\) and for D = 2r s , 

th F,s = ~ 2kD ' P D fa ' ln (i - Y F,s ) 
where D - diameter of evaporating droplet 


( 12 ) 


rearranging 


(13) 


(14) 
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Unity Lewis Number (Le) 

Assume Le = 1, it implies that the mass transfer rate is equal to the heat transfer rate, i.e. 

D FA =a g (15) 


where 


a 


g 


r k A 


pc 


p J 


(15a) 


where a g - thennal diffusivity of gas 
k - thennal conductivity 
c p - specific heat at constant pressure 


Then, 


Pg D fa ~ 



\ c p J 


Define the mass transfer number, B M '- 




y y 

^ F,drop ^ F,s 


Since Y Fi00 =0 and Y Fdrop = 1, equation (17) is simplified as: 


b m 


l F,s 


1 ~Y, 


F,s 


(16) 


(17) 


(18) 


Y F,s ~ 


1 + 


V P F,s 


-1 




n-1 


V M F J 


where Mp - molecular weight of fuel [kg/kg-mol] 

M a - molecular weight of air [kg/kg-mol] 

P - ambient pressure [kPa] 

Pp s - fuel vapor pressure at droplet surface [kPa] 

Rearranging In(l-Yp) in term of B\ t yields: 

ln(l -Y f )= -ln(l + B m ) 


(19) 


( 20 ) 
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Substituting (16) and (20) into (14) yields the rate of evaporation of a fuel drop at the surface: 


fh F,s 


■ 2nD 


' A 


\ c p J 


ln(l + B m ) 


( 21 ) 


Reference Conditions 

For better accuracy, the choice of values of k g and c Pig are evaluated at the following reference 
temperature and composition using the “one-third rule”: 


Since Y F oo ~ 0 , 


T =T +-(T -T ) 

r s 3 ' v ° 0 S 

(22) 

Y F,r ~ Yf,s + ^ _ Yf,s ) 

(23) 

7 -2y 

1 F,r ~ 

(24) 

Y — 1 -Y 

1 A,r 1 1 F,r 

(25) 


Fuel-air Mixture 

Therefore, the reference thermal conductivity and specific heat at constant pressure are estimated 
as: 


^=Y A , r -k A (T r )+Y Fr -k„,(T r ) (26) 

( ... = Y ,, )+ C pr , (T r ) (27) 

All properties for air and different fuels are provided in Appendix A and B. 


Evaporation Constant 

At steady-state period, the droplet diameter D at any instant may be related to its initial diameter 
Dq byD“-law: 


Dq —D 2 = A st t 


(28) 
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where A st [ ] is the steady state evaporation constant 

s 


4 = 




\ C p J 


811(1 + 5 ,,) 


Pf 


(28a) 


2 

The Z) -law states that the square of droplet diameter is a linear function of time where the 
evaporation constant apparently represents the slope of the equation. The larger the A st , the 
shorter the time it takes for the droplet to vaporize completely. 


Heat Transfer Number 

Consider conductive and convective heat fluxes across a thin shell surrounding the evaporating 
droplet, the heat transfer number is defined as the ratio of enthalpy available in the surrounding 
gas to the energy required to vaporize the fuel: 

D _ C P ,g(T~~T S ) pm 

°T ~ (t -T \ 29 

^ C p,dropV s * drop ) 


where L - latent heat of fuel vaporization corresponding to fuel surface temperature [ — ] 

kg 

For simplicity, one can neglect the energy required to raise the droplet temperature to the surface 
temperature. Then, equation (29) becomes: 


B T 


C p,g fco - T s ) 

L 


(30) 


When heat transfer dominates the evaporation process, the rate of evaporation of a fuel droplet at 
the surface is then described by: 


m F s = 2 nD 




\ c p ) 


ln(l + B t ) 


(31) 


Bt versus Bm 

Estimation of rate of fuel evaporation using equation (31) is only good for steady-state 
conditions. Nevertheless, equation (21) applies under all conditions, including the heat-up 
process of droplet (ref. 6). 
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However, under steady state conditions, B M = B r = B and equation (21) and (31) are identical. 
Therefore, droplet evaporation rate can be written as: 


m Fs = 2 ttD 




\ c p J 


ln(l + B) 


(32) 


2.2.2 Heat-up Process 

According to Chin 6 , serious error may be incurred in the calculation of fuel evaporation rate and 
droplet lifetime if the transient heat-up process is neglected. In fact, for many fuels at high 
ambient pressure and temperature, the transient heat-up process constitutes a significant portion 
of the droplet evaporation time. 

At the steady-state period, the heat used in vaporizing the fuel is given by: 


a 


= m F L 


Substituting equation (14) into (33) yields, 


Qss = 2 nD 


f k A 


\ c pj g 


ln(l + B m )L 


(33) 


(34) 


Including the heating process, the actual heat transfer is estimated as: 

Qac, = 2 XD( T . ~ T , ) s ln ^ ~ - 

B>, 


Then, the rate of change of the droplet surface temperature is given by: 

dT s _ Qact - Qss 
dt c pF m drop 

Substituting equation (34) and (35) into (36) and rearranging gives, 

dT s _ m F L ' B t | ' 

dt c p,F m drop \ d M ) 

where 

n n n 3 

m dmp ~~Pf d 


(35) 


(36) 


(37) 


(37a) 
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Note that 


d n j. 

m /• ~~r (~7 Pf d ) 
dt 6 


(38) 


Equating equation (14) and (38) and rearranging gives the rate of change of droplet size: 


dD 

dt 


41n(l + Sj 


( , A 


Pf D 


\ C PJ 


(39) 


2.2.3 Droplet Lifetime 

Since the rate of chemical reactions in many practical combustion systems are so high, the 
burning rate is mainly controlled by the fuel evaporation process. Therefore, droplet lifetime is 
important in such situations because it determines the residence time needed to ensure 
completion of combustion. 

Assume the final droplet diameter, D 0 , equal to zero and rearranging the D" -law, the steady state 
droplet lifetime is readily obtained by: 


l e,sl 


Do 

/L 


(40) 


2.2.4 Convective Effect 

All derivations shown above are only good for droplet at stationary condition. In addition, it is 
known that convection may enhance both mass and heat transfer during the evaporation process. 
Moreover, to include the convective effect into the evaporation rate equation is straightforward. 
Equation (23) is then modified by replacing the coefficient with the Nusselt number correlation 
(ref. 5): 


m F = NutjD 


r 


\ c p j 


l n (l + s„) 


where 


Nu = 2 + 0.6 Re°' 5 Pr°' 33 


(41) 


(41a) 


Re, 


pA-Ur'D 

Pa 


(41b) 
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All physical properties should be evaluated at the reference temperature, T r , except for p A , Pr A , 
and p A . Please refer to Appendix B for air properties. 


3. Numerical Method 

Six ordinary differential equations are to be solved for the six dependent variables x, z, Ud, w d , A 
and T s . The droplet trajectory is defined by the set of x and z values. A 4 th -Order Runge-Kutta 
explicit method 8 was used to solve these equations. The Runge-Kutta explicit method is an ideal 
numerical scheme for solving ordinary differential equations using Microsoft®" Excel software. It 
is a self-starting method with good stability characteristics. The time step-size can be changed as 
desired without any complications for higher-order schemes. 

There are totally six sets of coupled equations, namely the time rate change of x, z, u d , w d , D , and 
T s , along with their solutions, as shown below: (subscript n stands for the n' h time step) 


[1] 

= =/i(wrf) 

dt 

(1) 


x„+i =x n + Uk l +2k 1 + 2k i +k 4 ) 
6 

(la) 

where 

k \ =A t-f x (u d , n ) 



k 2 =A t-f x (u dn +y) 

k 3 =At-f l (u dn +y) 
k 4 = At • f x (u d n +/ 3 ) 

(lb) 

[2] 

^T f ^ w d =/ 2 ( w J 
dt 

(2) 


z n+ 1 = z n + — ( kz x +2kz 2 + 2 kz 3 + kz 4 ) 

(2a) 


6 


where 
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[3] 


where 


[4] 


where 


kz x = At ■ f 2 (w dn ) 
kz 2 = At • f 2 (w d n + 

kz 3 =At-f 2 (w dn +^~) 
kz 4 = At • f 2 (w d n + lz 3 ) 


(2b) 


du d 

dt 




)£/ 


"“gf? 






= A{ u d’ w d ’D,T s ) 


^d,n + 1 


« +t( / 1+2/2 +2^3 + ^) 


(3) 

(3a) 


A =^-/ 3 (u d> n,w d>n ,D n ,T s>n ) 

, * „ , /, /z, ^ md , _ mt i . 

^2 = At- f 3 (u dn + —,w d>n +—,D n +-^-,T sn + -^-) 

, * „ , U /z, ^ md 9 _ . 

^3 = At- f s {u dn +y,W rfj „ + -^,D n + 

h = At- / 3 (w rf , + / 3 , + fe 3 , D n + m d 3 , T + /h/ 3 ) 


(3b) 


dw d 

dt 


l P g ( w d ~ w g )pR \ A d C D + (/^ -PdVdS 


P d y d 


= M 


u d,w d ,D,T s ) 


(4) 


W 


df.W+l 


= A/,„ +^(/z,+2/z 2 +2/z 3 +/z 4 ) 
6 


(4a) 


lz\ =A t-f 4 (u dn ,w dn ,D n ,T sn ) 

, * /- / /| /z, ^ wJi 

/z 2 = At- f 4 (u dn + - , W rf> „ +—, D n + — 

, k r , h /z, ^ md 

lz 2 = A t-f 4 {u dn +j,w d ^ n +—,D n + — 

lz 4 = At- f 4 {u dn + l 3 , w dn + lz 3 , £>„ + mt/ 3 , 


- T -\ b 

’ s ’” 2 

L t +" 


+ /H/ 3 ) 


■) 

2 


) 


(4b) 
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[5] 


where 


dD _ A 
dt 2D 


fs(D,T s ) 


D n+l = D n + —{md x +2md 2 + 2md 3 + md 4 ) 
6 


md x = A t-f 5 (D n ,T sn ) 

, . ^ „ md, _ mt, „ 

md 2 = At ■ f s (D n + — , T s >fl + — ) 

, . „ . ^ md-y mt 2 . 

/«J 3 = At- f 5 (D„ + ~y~ , T s ,n + -^) 

md 4 = At ■ f 5 ( D n + md 3 , T s n + /h/ 3 ) 


(5) 

(5a) 

(5b) 


[ 6 ] 


where 


dT s _ m F L ' B r | 

^ C p,F m drop , 




T —T 

s,n + 1 


l 

+ — 
6 


(mt l +2mt 2 +2 mt 3 +mt 4 ) 


( 6 ) 

(6a) 


mt x = At -f 6 (u dn , w d n , D n , T s n ) 

„ , 6 /z, znJ. _ /«/, _ 

mt 2 = A t- f 6 (u d n + —,w d n + — ,D n + -^-,T s n +-^-) 

. . , l 2 lz 2 md 2 _ x 

mt 3 =At-f 6 (u dn +^-’ w d , n + ^’ D » + ~2T Js - n + ~f ) 

mt 4 = At ■ f 6 (i u d n + / 3 , + /z 3 , £>„ + md 3 , r s>JI + mf 3 ) 


Only every n th cycle (as specified by the user) is saved for plotting. This greatly saves on storage 
and increases the speed of post processing. We have chosen to enter the data in SI units in the 
unlocked cells. The required conversions are done in the locked cells. When the user becomes 
familiar with the spreadsheet, the spreadsheet can be unlocked as there is no password and the 
user can adapt the spreadsheet as required. 
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The report can be accessed on the web at: 

http://gltrs.grc.nasa.gov/reports/2004/TM-2004-2 129 1 0/TM-2004-2 129 1 0.html 

along with the interactive Microsoft® Excel spreadsheet 1 that computes and plots 

data with and without spray evaporation can be accessed through hyperlinks located on the back 

of the title page and on page 17, or directly on the web at: 

http://gltrs.grc.nasa.gov/reports/2004/TM-2004-212910/SprayEvapVBONOFFSept23.xls 


The interactive Microsoft® 1 Excel spreadsheet is also available on CD-ROM as a separate 
document. Additional copies of the Microsoft® 1 Excel spreadsheet can be requested by e-mailing: 
Dan.L.Bulzan@nasa.gov . 

The CD-ROM also contains NASA/TM — 2002-211710 and supplemental Microsoft® Excel 
interactive spreadsheet that computes and plots data without spray evaporation. 


'To access the interactive spreadsheet, Microsoft 8 Excel 2002 or higher is recommended to view the fde and for 
proper functionality. It can be opened through your browser, however, saving to the hard drive is recommended. If 
you cannot access this fde, please contact: Dan.L.Bulzan@nasa.gov . 
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4. Equations Summary for Tracking an Evaporating Spray in a Crossflow 


Droplets 

Properties: 

Initial conditions: 

Crossflow 

Properties: 

Ambient conditions: 
Numerical method. 

Ballistics: 

Evaporation: 


Inputs 


P~32> Ad, V d, Pd.288.25Kt T cr it, T } m , Cexp, Lpb n , IVl W /, a, b 
Xo, Udo, Z 0 , Wdo, T init , 


Mg, U g , W g 
Tec, P, g 

At, n (total cycle number), iidata, £ 


Outputs 

x, Ud, z, Wdo, Cdt 

D, T sur f 
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Equations 


Mass Transfer Number: 


If (Group 1 fuel - ) then 


P F,s = P vap X ) = exp 


' b ' 

a — - 


V 7,-43, 


[kPa] 


else 


End if 


l°gio P vap ~ A + — + C ^°^ 10 Ts + D p s + pp s 

* S 


P F,s = P va P X ) [nunHg] 


y FJ ,= 


1+ 


- — i 

p„ 

V F ’ s J 


V M, A 


-|-1 


F J 


B 


Y, 


F.s 


M 


1 -Y, 


F.s 


Reference Conditions: 


K =r,+t(7’.-r,)[K] 


v = — Y 

1 F.r l ) lF ’ S 


Y A =l-Y Ftr 


Thermal Properties of Air 


k a =2-10 


' n T 3 -5-10 


-8 rri 2 


Tf+ 0.0969-10 


~ 3 T r +0.8289 1CT 3 [ 


W ■ 
m ■ K ' 


“Group 1 fuels - DF 2, JP 4, JP 5 and n-Heptane 
Group 2 fuels - Jet-A (C 12 H 23 ) and Water (H 2 0) (please refer to Appendix A for more details.) 
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(10) 


c pA =-5 10 ' n r ; 3 + 2 -\ 0 ~ 7 T t 2 -llO~ 5 T r + 1 . 004 1 [— ] 

kg 

fl A =( 710 " 8 ri- 0 . 0003 ri + 0 . 6227 r TO + 18 . 761 )- 10 ' 7 [^] ( 11 ) 

m 1 


Pr A = 2 - 10 “ 13 ri - 10 “ 9 ri + 2 ■ lO -6 ^ - 0 . 00097 L + 0.8632 ( 12 ) 


Pa = 355 . 91 rj L0032 A ( 13 ) 

m J 


Thermal properties of hydrocarbon fuel: 

If (Group 1 fuel) then 


Else 


Pf ~ Pf, 288 . 6 K 


1 - 1 . 8 C CT (r- 288 . 6 )- 0.09 ^ 288 ' 6 ^ 


(T crit - 288 . 6 )" 


m 


/ \ kJ 

c p .r, = ( 0.363 + 0 . 000467 r r )( 5 - 0 . 001 /), >2 „. 6i .)[— ] 

kgK 


c p ,f = ( 0.76 + 0 . 00335 r,„,, ((O.OOl/v f 5 [-^] 


« = 2 - 0.0372 


F T \ 2 

1 r 

\ T bnJ 


k Fv = 10 “ 6 [l 3 . 2 - 0 . 0313 ( 7 ),„ - 273 )] 


Ft \ n 


V 273 y 


w 
[ — ] 
m-K 


L = L 


Tbn 


T ■ -T 

cnt s 

T T 

V 1 crit 1 bn J 


\ 0.38 




- 1 


L drop j 

Pf=AB v Tc ' [-^j] 
cm 


( 14 ) 

( 15 ) 

( 16 ) 

( 17 ) 

( 18 ) 

( 19 ) 

( 20 ) 
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( 21 ) 


k Fv =A + BT r + CT r 2 [-^—] 

m ■ K 

C„ Fv =A + BT r + CT; + DT? + ET } 4 [ - ] 

p ■ mol-K 

Cp,F = A + BT drop + CTdrop + DT drop [ /;;o / _ g ^ 


L = A\ 


Tdrop yl 


[-^] 

mol 


End if 

Thermal properties of fuel-air mixture: 


k s =Y A ,-k A {T r )+Y F ,-k v {T r ) 


C p,g = Y A, r ■ C p,A ( T r ) + Y F,r ' C p,v i T r ) 


Droplet trajectories: 


Ur = J(" d -U g ) 2 +(w rf -W g ) 2 


Re rf =■ 


2 P, 


U t 


d.s 


M, 


C D =< 


24 


Re, 


, , 1 „ 2/3 

l + -Re d 
6 


0.424 


Re rf <1000 


Re rf > 1000 


P^^r = ~P g (« d -« g ) P 


dt 2 


g\“d ”g) R 


djC D 


pF J ^r=~p g («’ J -wj \uU d c D +(p- Pj )v d g 


dt 2 


122 ) 

(23) 

(24) 


(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 
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( 32 ) 


Convective effect: 


Heat transfer number: 


Evaporation constant: 


Heat-up process: 


Evaporation time: 


Nu =2 + 0.6 Re' 1 / Pr*', 33 


in F = NujiD — ln(l + B m ) 


B t = 


C p „ (71 - T s ) 


k_\ 8/n(l + B m ) m 2 


T/’ff’lkg] 


dT s _ m i L ( B_ T J ^ 

dt C p,F m drop { B M ) S 


dD A m 

= [-] 

dt 2D s 


D 2 


(33) 

(34) 

(35) 

(36) 

(37) 

(38) 

(39) 

(40) 

(41) 
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5. User’s Manual 


This program is written in Microsoft® Visual Basic Excel. There are six sheets in the program, 
namely the ‘CoverPage’ sheet, ‘Instruction’ sheet, ‘Process’ sheet, ‘Data’ sheet, ‘Trajectory’ 
sheet, and ‘Evap’Code sheet. 

‘CoverPage’ sheet 

Relevant information about authors is provided in this sheet. If there are any questions or 
comments regarding this program, please feel free to contact us. 

‘Instruction’ sheet 

The instruction sheet contains a description of each sheet in the program, as well as the user’s 
manual. 



Figure 5.1 - ‘Process’ sheet 
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‘Process’ sheet 

The process sheet contains all user-inputs necessary for performing computations. First of all, the 
user needs to choose a fuel type from the pull down menu. This program contains built-in a 
database of properties of each listed fuel necessary for evaporation computations (see 
Appendix A - Fuel Properties). Other than that, the user can specify any hydrocarbon fuel not 
listed here, provided that the properties of each specified fuel are available. 

The cells highlighted in light-brown are the user inputs. The cells highlighted in pink contain 
computed values associated with the light-brown cells; therefore, they are locked to prevent the 
user from modifying their values. When values are entered into the fonnula cells, the formulas 
are erased and linkages to other cells are interrupted. That is why for this version we chose to 
lock the closed cells without a password. 

One feature added to the previous work (ref. 1) is the ‘Evaporation ON/OFF’ switch. If the user 
chooses to turn off the ‘Evaporation switch’, the Microsoft'® Excel program will run as if there is 
no evaporation and will give the solution for the droplet ballistics only. 

If the ‘Evaporation switch’ is turned on, one main concern is the computational time because it is 
considered very computation-intensive to use Runge-Kutta method in a Microsoft® Excel 
spreadsheet for tracking a single evaporating droplet in crossflow. The computational time varies 
mostly with the CPU power available and the total cycle number as well as the time step size. To 
remedy this problem, four options have been proposed and it is up to the user to select one of the 
four options before running the code: 

a. Option 1 - Properties are calculated four times each cycle (more accurate but very slow) 

b. Option 2 - Properties are calculated one time each cycle (less accurate but fast) 

c. Option 3 - Properties are calculated four times during droplet heating. After droplet 
reaches steady state temperature, properties are not calculated any more. The values are 
used from the last cycle of droplet heating. This option requires user to input a constant 
value of epsilon, e. It defines the condition when the steady state temperature will occur. 
(dT/dt = e) 

d. Evap time estimate - Before computation, it’s always good to estimate the steady state 
evaporation time (droplet lifetime) of the droplet of size ryj- Click the Update button to 
start the estimation. The computational speed depends mostly on the droplet size r 32 and 
ambient conditions. For example, a Jet- A droplet of 50-micron diameter (r 32 =25 micron) 
is estimated to completely vaporize within 0.004s. Using At=1.0e-7s, at least n=40,000 
cycles are required. 

4 th -Order Runge-Kutta method calculates temperature derivative (dT/dt) four times in each cycle. 
Option 3 will force the program to check all four temperature derivatives in each cycle. If one of 
the temperature derivatives in a cycle is less than the constant epsilon, properties calculation will 
not be performed afterward. Based on the experiment, epsilon between 0.1 and 1 is good enough. 
For instance, a case has been perfonned using two constant epsilons, say 0.0001 and 0.1. 
Solutions from both epsilons yield similar results. Note that this may not apply to cases with fuel 
with properties highly sensitive to temperature. 
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From figure 4, the last two user inputs, Total Cycles (C32) and Data taken every ## cycles ( C33 ) 
are the important features that were added for monitoring the amount of output data. The 
number assigned in pink cell “E34”must be kept below 65,536; the cell will turn red if this 
condition is not satisfied. Keeping the value below the limit can be done by changing the value 
in the cell “C35” 

After all the inputs have been specified, clicking the Update button will instruct the program to 
update the data in the ‘Data’ sheet as well as all plots in the ‘Trajectory’ and ‘Evap’ sheets. 

In summary, the user needs to do the following steps to run the program: 

1. Go to the ‘Process’ Sheet click on the Process tab. 

2. Enter input values in the light-brown cells. 

3. Adjust the value in the cell C33 so that the computed value in cell E34 is less than 
65,536. 

4. Choose the fuel type from the ‘ Fuel ’ option menu. 

5. Choose the ON/OFF evaporation switch from the evaporation option menu. 

6. If Evaporation is ON, choose one of the four options: Option 1, 2, 3 or Evap time 
estimate. 

7. If ‘ Evap time estimate ’ is selected in step 6, click the Update button. The code will 
estimate evaporation constant and time. 

8. Based on the estimated value of t stevap (droplet lifetime), adjust and input appropriate 
’Time step At' and ’Total cycles n'. 

9. Click the Update button. 

10. Observe the droplet profiles on the solution plots at both ‘Trajectory’ and ‘Evap’ 
sheets. 

11. Repeat step 1 through step 10 for different input values 

12. Click the Clear Data button to clear the data (Optional). 

An additional feature in this program is the option to store the computed data into a TecPlot 
format file. This feature provides the user a flexibility to plot the data using other software such 
as TecPlot. In summary, user needs to do the following steps to save the data into TecPlot format 
file: 

1. Select "Yes" to save the data into a file 

2. Specify the path and the filename to save into a file (e.g. C:\Document\result.dat) 

3. Select the unit length (i.e. meter or millimeter) 

4. Click the Update button 

Other than computational time, another concern is the computer memory usage because the code 
needs a lot of memory space for those six variables arrays to be solved. As discussed above, the 
feature of estimating the steady state droplet lifetime enables us to predict the total cycle number 
n required to complete a whole evaporation process for each fuel droplet. By making use of this 
advantage, a so-called ‘dynamic array size ’ is added to the existing code such that the allocated 
memory space for each computation will closely follow the total cycle number n entered by user. 
Using the example in ‘ Evap Time Estimate’’ option where n=40,000, each variable array size will 
be set to «+2=40002. 
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‘Data’ sheet 

This sheet contains the solution data computed by the program. All the solution data will be 
listed separately at eight columns A to H for plotting. Table 1 shows the corresponding variables 
of each column: 


Column 

Variable 

Units 

Definition 

A 

Time 

sec 

Time 

B 

X 

m 

Droplet trajectory in x-direction 

C 

Z 

m 

Droplet trajectory in z-direction 

D 

u d 

m 

V 

Droplet velocity in x-direction 

E 

w d 

m 

— 

s 

Droplet velocity in z-direction 

F 

C D 

n/a 

Drag coefficient 

G 

D/D n 

jim 

Nonnalized droplet diameter square 

H 

T s 

K 

Droplet surface temperature 


Table 5.1 - Variables with their corresponding column in ‘Data’ sheet 


‘Trajectory’ sheet 


As shown in figure 5, there are five graphs on this sheet, namely droplet trajectory, droplet 
velocity profile, drag coefficient, Cd, as a function of time, droplet velocity profiles as a function 
of time, and droplet trajectory profiles as a function of time. 


NASA/TM— 2004-212910 


26 






[2 Microsoft Excel - Spray_EvapVB0N0FFSept21 




File Edit View Insert Format Tools Data Window Help 


J5|xJ 

Type a question for help ..fix 



Ready 

1 Start! 


NUM 


] Microsoft Excel - Spra... 


U 9 11:09 AM 


Figure 5.2 - ‘Trajectory’ sheet 
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‘Evap’ sheet 


There are two graphs on this sheet, as shown in figure 6, namely droplet size and surface 
temperature histories. It also shows the values of the estimated steady-state evaporation constant 
and evaporation time. 


E Microsoft Excel - Spray_EvapVB0N0FFSept21 


f) File Edit View Insert Format lools Data Window Help Type a question for help . _ 5 x 



* Start] | [x] Microsoft Excel - Spra... 11:13 AM 


Figure 5.3 - ‘Evap’ sheet 
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6. Results and Discussions 


To understand the performance of this Excel code, several results have been obtained and 
displayed in the following pages using the input parameters as shown in Table 6.1. Only the 
trajectory and the time variation of normalized droplet size and surface temperature are selected 
as the results to be discussed. 


P=1 atm.T=500K.R32=25|im 



6.1(a) Normalized diameter square vs time 

P=1 atm.T=500K.R32=25|xm 


500 
480 
460 
440 
420 
' 400 
380 
360 
340 
320 
300 


Water 

DF 2 

Jet- A 


0.1 


0 0.02 0.04 0.06 0.08 

Time (s) 

6.1(b) Droplet surface temperature vs time 

P=1atm. T=500K. R32=25mm 


0.12 0.14 



6.1(c) Trajectory of fuel droplet after injection 


P=10atm,T=500K,R32=25[im 



6.2(a) Normalized diameter square vs time 


P=10atm,T=500K,R32=25|im 



6.2(b) Droplet surface temperature vs time 


P=10atm. T=500K. R32=25mm 



6.1(c) Trajectory of fuel droplet after injection 


Figure 6.1- Injection of fuel droplet at 
T=500K and P=latm 


Figure 6.2 - Injection of fuel droplet at 
T=500K and P=10atm 
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Ts (K) D A 2 D0 A 2 


P=1 attn.T =1000K.R32=25jim 



6.3(a) Normalized diameter square vs time 
P=1atm,T=1000K,R32=25|im 



Time (s) 

6.3(b) Droplet surface temperature vs time 


P=10atm,T=1000K,R32=25|im 



6.4(a) Normalized diameter square vs time 



Time (s) 

6.4(b) Droplet surface temperature vs time 


P=1«m. T=1000K. R32=25mm P=10atm, T=1000K, R32=25mm 




6.3(c) Trajectory of fuel droplet after injection 6.4(c) Trajectory of fuel droplet after injection 


Figure 6.3 - Injection of fuel droplet at 
T=500K and P=latm 


Figure 6.4 - Injection of fuel droplet at 
T=500K and P=10atm 
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In figure 6.1(a), water droplet of diameter 50 micron at 500K 1 atm takes a much longer time to 
reach complete vaporization i.e. around 0.13s than Jet-A and DF 2 do. As shown in figure 6.2(a), 
increasing the ambient pressure to 10 atm further extends its droplet lifetime to around 0.16s, as 
known from the general physics results. Similar conclusions will be obtained when comparing 
figure 6.3(a) with 6.4(a) except that each fuel droplet lifetime is much shorter at such a high 
ambient temperature (T=1000K). 

Figures 6.1(b) and 6.2(b) show the time variation of droplet surface temperature at the same 
ambient temperature i.e. 500K and two different ambient pressures, i.e. 1 atm and 10 atm 
respectively. In figure 6.1(b), each droplet reaches its steady state temperature and it never goes 
beyond its normal boiling point. While at high pressure as shown in figure 6.2(b), the boiling 
point increases with pressure and consequently the steady state temperature is higher. These 
observations also apply to figures 6.3(b) and 6.4(b). 

The droplet trajectory after injection and before complete vaporization for each droplet is shown 
in figures 6.1(c), 6.2(c), 6.3(c) and 6.4(c). Due to the decreasing droplet mass (as a result of 
vaporization) followed by the momentum loss, the downward crossflow forces the evaporating 
droplet to drop almost vertically (or in z-direction) at the very end of its drop life. Water, with a 
higher heat capacity, therefore penetrates deepest across the crossflow among all the fuels. 


Input 

Value (unit) 

Input 

Value (unit) 

t*32 

25 jam 

w g 

-38 m/s 

X 0 

0m 

To. 

1000 K 

c 

a 

o 

-2.4 m/s 

P 

lOIOkPa 

Zo 

0m 

g 

Om/s 2 

Wd Q 

0 m/s 

At 

1 .00E-07 sec 

T init 

300 K 

n 

100000 Cycles 

M a 

28.97 kg/kmol 

8 

1 K 

U, 

0 m/s 




Table 6.1- Input parameters used in discussion 
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Appendix A — Fuel Properties 


In this code, two different groups of fuels are considered. The reason is because the thermal 
properties for these two groups of fuels are referring to different sources as shown below. 


Group 

Source 

Fuels 

1 

Chin (refs. 3 and 6) 

DF 2 (diesel oil) 
JP 4 
JP 5 

n-Heptane 

2 

Yaws (ref. 7) 

Jet-A (C 12 H 23 ) 
Water (H 2 0) 


Table A. 1 - two groups of fuels 


Fuel/properties 

DF 2 

Jet-A 

(Ci 2 H 23 ) 

JP 4 

JP 5 

n-Heptane 

Water 

(H 2 0) 

Pd, 288 . 6 K (kg/m 3 ) 

846 

915.91 

773 

827 

687.8 

1036 

Tcrit(K) 

725.9 

737.0 

612 

648.8 

540.17 

647.3 

T b „ (K) 

536.4 

529.0 

420 

495.3 

371.4 

373.16 

C«d(1/K) 

0.00046 

n/a 

0.000557 

0.000485 

0.000715 

n/a 

Ljbn (kJ/kg) 

254 

341.200 

292 

266.5 

317.8 

2191 

M d (kg/kmol) 

198 

181.321 

125 

169 

100.16 

18.02 

a* 

15.5274 

n/a 

15.2323 

15.16 

14.2146 

n/a 

b* 

5383.59 

n/a 

3999.66 

4768.77 

3151.68 

n/a 


*Clasius-Clapeyron constants 


Table A. 2 - Physical properties of fuels (Group 1 and 2) 

All fuel properties are calculated in subroutines transfernumber(Ts) and 
fuel_thermal_prop (Ts) where Ts refers to surface temperature of a fuel droplet. 
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Group 1 Fuels: DF 2 (diesel oil), JP 4, JP 5, and n-Heptane 


At any given droplet surface temperature, T s , fuel vapor pressure at droplet surface, P FyS , is 
approximated by Clasius-Clapeyron equation (or called Antoine equation (ref. 4)) [kPa]: 


p i ,s = P va P ( r s ) = CX P 


b 


T - 43 


(A.l) 


where a and b are constants for certain fuels (provided in Table A. 2 shown above). 
T s - droplet surface temperature [K] 


Meanwhile, variation of k Fv 
by Chin (ref. 3): 


W 


m ■ K 


] and c 


p,Fv 


— ] of fuel vapors with temperature are given 
kg 


k Fv = 10“ 6 [l3.2- 0.03 l3{T bl , -273)] 


f T \ H 

1 r 

V273y 


(A.2) 


c p , Fv = (0-363 + 0.0004677) )(5 - 0.001^ >28JL6JC ) 


(A. 3) 


where 


n = 2 — 0.0372 


r T \ 2 

1 r 

\ T bn J 


(A. 3a) 


Pf 288 . 6 k " density in kg/m 3 at 288. 6K. 


3 

Liquid fuel density [ — - ] is needed to calculate A st : 


m 


Pf=P 


F,2&&.6K 


(A.4) 


where C exp is the expansion coefficient (provided in Table A.2). 
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Latent heat of fuel vaporization [ — ] corresponding to fuel surface temperature is estimated by: 

kg 


L ^ Tbn 


( 



T -T 

V crit bn J 


(A. 5) 


Ljbn ~ Latent heat of fuel vaporization at normal boiling point [ — ] 

kg 


Specific heat at constant pressure of fuel liquid is given by [ ]: 

kgK 


c„. r = (0.76 + 0.003357;„„)(0.001p,) 05 


(A.6) 
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Group 2 Fuels: Jet-A (C12H23) and Water (H 2 0 ) 

All equations for the fuel properties can be found in Yaws (ref. 7), as shown below: 


Fuel liquid density [ 


g 

3 

cm 


]: 


Pf — AB 



(A. 7) 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

0.29292 

0.3471 

B 

0.26661 

0.274 

n 

0.298 

0.2857 


O’ 

Table A.3 - Coefficients for Fuel Liquid Density [ , ] in Group 2 fuel. 

cni 


Latent heat of vaporization of fuel liquid [ 


kJ_ 

mol 


]: 


r 

L = A 1 

v 



A 


J 


n 


(A. 8) 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

73.509 

52.023 

n 

0.347 

0.321 


Table A.4 - Coefficients for Latent Heat of Vaporization of fuel liquid [ ] in Group 2 fuel. 

mol 
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Specific heat at constant pressure of fuel liquid [ ]: 

mol — K 


C p F — A + B + CT, 


drop 


[ drop 


+ DT 


drop 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

142.238 

92.053 

B 

1.5261 

-0.039953 

C 

-0.0034477 

-0.00021103 

D 

0.0000032968 

0.0000005347 


Table A. 5 - Coefficients for Specific Heat at constan t pressure of fuel liquid [ ] 

mol — K 

Group 2 fuel. 


Specific heat at constant pressure of fuel vapor [ ]: 

mol — K 

C p Fv - A + BT r + CT; + DT , 3 + ET? 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

-128.032 

33.933 

B 

1.4622 

-0.0084186 

C 

-0.00086193 

0.000029906 

D 

0.00000018462 

-0.000000017825 

E 

0.00000000000036227 

3.6934E-12 


Table A. 6 - Coefficients for Specific Heat at constant pressure of fuel vapor [ ] 

mol — K 

Group 2 fuel. 


Thermal conductivity of fuel vapor [ 


W 


m ■ K 


]: 


kp v — A + BT r + CT r 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

-0.01184 

0.00053 

B 

0.000061839 

0.000047093 

C 

0.000000025082 

0.000000049551 


Table A.7 - Coefficients for Thermal Conductivity for fuel vapor [ ] in Group 2 fuel. 

m ■ K 


(A. 9) 


in 


(A. 10) 


in 


(A. 11) 
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Vapor pressure [mmHg]: 


logio P vap =^ + ^r + C log 10 T s + DT S + ET ; 


(A. 12) 


Then, let 


F,s 



(A. 12a) 


Coefficients 

Jet-A (C 12 H 23 ) 

Water (H 2 0) 

A 

-50.5512 

29.8605 

B 

-2705.3 

-3152.2 

C 

28.273 

-7.3037 

D 

-0.045702 

0.00000024247 

E 

0.000020443 

0.000001809 


Table A.8 - Coefficients for Vapor Pressure [mmHg] in Group 2 fuel. 
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Appendix B — Air Properties 


Variation of air properties with temperature (range 100K through 2000K) can be obtained from 
correlation provided in Incorpera 5 . Air properties are all calculated in subroutine air () . 

W 

Thermal conductivity [ ] : 

in ■ K 


k A = 2-10 ~ U T , 3 - 5 10 ~*T; - 0.0969 10“ 3 r,. + 0.8289 10“ 3 

kj 

Specific heat at constant pressure [ — ]: 

kg 

c p , A =-5-io -11 r 3 + 2-io~ 7 r r 2 -i-io“ 5 r r + 1.0041 


(B.l) 


(B-2) 


Viscosity [ — 
m 2 

n A =(7io _8 ri -o.ooo3ri +0.6227^ + 18.761)- io -7 

Prandtl number: 


Pr 4 = 2 ■ 10“ 13 rj -10“ 9 ri + 2 ■ 10 ~ 6 T 2 - 0 . 00097 L + 0.8632 


ks 

Density [— y ]: 

m 3 


(B-3) 


(B-4) 


Pa 


355.91 T 


- 1.0032 


(B-5) 
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Appendix C — Microsoft® Visual Basic Code 

In the Microsoft® Excel spreadsheet, select ‘Tools->Macro->Visual Basic Editor’, a new 
window will open showing all the numerical codes in different sheets or modules. 

Main Code 


Private Sub CommandButton3_Click() 

'Declare ballistics variables 

Dim n As Long, nn As Double, nnnn As Long, speedup As Long 
Dim Cdm As Double, Rems As Double 
Dim dtt As Double 

Dim kl As Double, k2 As Double, k3 As Double, k4 As Double, 11 As Double, 12 As Double, 13 As Double, 14 
As Double 

Dim kzl As Double, kz2 As Double, kz3 As Double, kz4 As Double, lzl As Double, lz2 As Double, lz3 As 
Double, lz4 As Double 
Dim Urm As Double 
Dim a As Double, b As Double 
Dim location 
Dim sf 

Dim at, bl, cl, dl, el, ft, gl, hi 
Dim unitconv 

'Adding these inputs for drop evaporation '(3/26/03) 

Dim mdl As Double, md2 As Double, md3 As Double, md4 As Double 
Dim mtl As Double, mt2 As Double, mt3 As Double, mt4 As Double 
Dim evap As Double 

'Variable time 

Dim begin, last, switch 

begin = Second(Time) + Minute(Time) * 60 + Hour(Time) * 3600 
switch = 0 


'Clear data 
Call Macro2 
Halt = False 

Msg = "Do you want to continue ?" 

Style = vbYesNo 

Style 1 = vbOKOnly 

Title = "Jet Flow in CrossFlow" 

Titlel = "Evaporation Time Estimate" 

Msg3 = "Drop vaporizes completely!" 

Msg size = "Please reduce droplet size!" 

Msg sizel = "Can't estimate. Drop size too big!" 

Call taperd 

'Dynamic array size 
asize = kl + 2 

ReDim xx(0 To asize) As Double, xp(0 To asize) As Double 
ReDim zz(0 To asize) As Double, zp(0 To asize) As Double 
ReDim D(0 To asize) As Double, Ts(0 To asize) As Double 
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'Sheets("Process").Cells(9, 10) = kl 
'Sheets("Process").Cells(10, 10) =asize 
If flag = 4 And r > 0.0001 Then 

Response = MsgBox(Msg_size, Style 1, Title 1) 

Cells(5, 3). Select 
GoTo 50 
End If 

Pi = 22# / 7# 

Worksheets("Process").CommandButtonl. Width = 0 
Worksheets("Process").CommandButtonl. Visible = True 

n = 0 'used for array value index 
nn = -1 'used for writing data into cell 
nnn = 2 'used for writing into sequence cell 
dtt = 0 

evap = 0# 

nl = 0 'used to keep track the time level recorded when drop vaporizes completely 

ss = 0 'used to turn on/off switch if temperature has reached steady 

sst = 0 'used to print steady state time 

speedup = Intfkl / 50) 'used to update the percentage bar 

'initial position and velocity 

xx(n) = xxi 

xp(n) = xpi 

zz(n) = zzi 

zp(n) = zpi 

'Assume uniform temperature within the drop 
Ts(n) = Tinit 

D(n) = 2# * r 
D02 = D(n) A 2 

If flag = 3 Or flag = 4 Then Sheets("Process").Cells(18, 10) = "Temperature is not steady" 

'Determine fuel type then select one of two approach to calculate fuel properties 
' 0 - Lebefvre's, 1 - Yaws' 

fueltype = Sheets("Process").ComboBoxl. Value 
Select Case fuel type 
Case "DF 2" 
otherfuel = 0 
Case "Jet-A (C12H23)" 
otherfuel = 1 
Apv = -50.5512 
Bpv = -2705.3 
Cpv = 28.273 
Dpv = -0.045702 
Epv = 0.000020443 

t 

Arho = 0.29292 
Brho = 0.26661 
nrho = 0.298 
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Akfv = -0.01 184 
Bkfv = 0.000061839 

Ckfv = 0.000000025082 
! 

Acpfv = -128.032 
Bcpfv = 1.4622 
Ccpfv = -0.00086193 
Dcpfv = 0.00000018462 

Ecpfv = 3.6227E-13 
! 

Acpf = 142.238 
Bcpf = 1.5261 
Ccpf= -0.0034477 

Dcpf= 0.0000032968 
» 

Alat = 73.509 
nlat = 0.347 

Case "JP 4" 
otherfuel = 0 
Case "JP 5" 
otherfuel = 0 
Case "n-Heptane" 
otherfuel = 0 
Case "Water" 
otherfuel = 1 
Apv = 29.8605 
Bpv = -3152.2 
Cpv = -7.3037 
Dpv = 0.00000024247 

Epv = 0.000001809 
» 

Arho = 0.3471 
Brho = 0.274 
nrho = 2# / 7# 

t 

Akfv = 0.00053 
Bkfv = 0.000047093 

Ckfv = 0.000000049551 
» 

Acpfv = 33.933 
Bcpfv = -0.0084186 
Ccpfv = 0.000029906 
Dcpfv = -0.000000017825 

Ecpfv = 3.6934E-12 
! 

Acpf =92.053 
Bcpf =-0.039953 
Ccpf= -0.00021103 

Dcpf= 0.0000005347 
! 

Alat = 52.023 
nlat = 0.321 
Case Else 
otherfuel = 0 
End Select 
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If onoff flag = 0 Then 
Tref = Tinit + (Tinf - Tinit) / 3# 

'air: viscosity 

mu = 0.00000007 * Tref A 3 - 0.0003 * Tref A 2 + 0.6227 * Tref + 18.761 
mu = mu * 0.0000001 'Ns/m A 2 

If otherfuel = 0 Then 
'fuel: liquid density 

rhof = rhofr * (1# - 1.8 * cexp * (Tinit - 288.6) - 0.09 * ((Tinit - 288.6) / (Tcrit - 288.6)) A 2#) 'kg/m A 3 
Else 

'fuel: liquid density 
Tdocrit = Tinit / Tcrit 

rhof = 1000# * Arho * Brho A ((-1#) * (1# - Tdocrit) A nrho) 'convert g/cm A 3' to 'kg/m A 3 
End If 

'air: density 

rhog = 355.91 * Tref A (-1.0032) 'kg/m A 3 
End If 


Sheets("Process").Cells(14, 12) = "Calculating" 

20 

dtt = dtt + dt 

If flag = 1 Or flag = 2 And onoff_flag = 1 Then Call driver_for_runge(Ts(n)) 

If flag = 3 And ss = 0 And onoff_flag = 1 Then Call driver_for_runge(Ts(n)) 

If flag = 4 And ss = 0 And onoff flag = 1 Then Call driver_for_runge(Ts(n)) 

Call Runge(xp(n), zp(n), D(n), Ts(n)) 

If (SuperExit) Then GoTo 50 

kl = dt * xp(n) 'replaced xp(n) 3/26/03 

11 = dt * A_x 

kzl = dt * zp(n) 'replaced zp(n) 3/26/03 
lzl = dt * A_z 
If onoff_flag = 1 Then 
mdl = dt * (-0.5) * lambda / D(n) '3/28/03 
mtl = dt * ct 

If ((D(n) + 0.5 * mdl) < 0.00000001) Then 
evap =1# 
nl = n 
GoTo 15 
End If 
End If 

If flag = 1 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + 0.5 * mtl) 

If flag = 3 And ss = 0 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + 0.5 * mtl) 
Call Runge(xpfn) + 11 / 2#, zp(n) + lzl / 2#, D(n) + 0.5 * mdl, Ts(n) + 0.5 * mtl) 

If (SuperExit) Then GoTo 50 

k2 = dt * (xp(n) + 11 / 2#) 'replaced xp(n) 3/26/03 

12 = dt * A_x 

kz2 = dt * (zp(n) + lzl / 2#) 'replaced zp(n) 3/26/03 
lz2 = dt * A z 
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If onoffflag = 1 Then 

md2 = dt * (-0.5) * lambda / (D(n) + 0.5 * mdl) '3/28/03 
mt2 = dt * ct 

If ((D(n) + 0.5 * md2) < 0.00000001) Then 
evap =1# 
nl = n 
GoTo 15 
End If 
End If 

If flag = 1 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + 0.5 * mt2) 

If flag = 3 And ss = 0 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + 0.5 * mt2) 
Call Runge(xp(n) + 12 / 2#, zp(n) + lz2 / 2#, D(n) + 0.5 * md2, Ts(n) + 0.5 * mt2) 

If (SuperExit) Then GoTo 50 

k3 = dt * (xp(n) + 12 / 2#) 'replaced xp(n) 3/26/03 

13 = dt * A x 

kz3 = dt * (zp(n) + lz2 / 2#) 'replaced zp(n) 3/26/03 
lz3 = dt * A_z 
If onoff_flag = 1 Then 

md3 = dt * (-0.5) * lambda / (D(n) + 0.5 * md2) '3/28/03 
mt3 = dt * ct 

If ((D(n) + md3) < 0.00000001) Then 
evap =1# 
nl = n 
GoTo 15 
End If 
End If 

If flag = 1 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + mt3) 

If flag = 3 And ss = 0 And onoff_flag = 1 Then Call driver_for_runge(Ts(n) + mt3) 

Call Runge(xp(n) + 13, zp(n) + lz3, D(n) + md3, Ts(n) + mt3) 


If (SuperExit) Then GoTo 50 

k4 = dt * (xp(n) + 13) 'replaced xp(n) 3/26/03 

14 = dt * A x 

kz4 = dt * (zp(n) + lz3) 'replaced zp(n) 3/26/03 
lz4 = dt * A_z 
If onoff flag = 1 Then 

md4 = dt * (-0.5) * lambda / (D(n) + md3) '3/28/03 
mt4 = dt * ct 
End If 

If flag = 3 Or flag = 4 And ss = 0 And onoff_flag = 1 Then 
If Abs(ct) < sstemp Then ss = 1 
End If 

If flag = 3 And ss = 1 And sst = 0 And onoff_flag = 1 Then 
Sheets("Process").Cells(18, 10) = "Temperature is steady at t = " & dtt * 1000# & "milliseconds" 
sst = 1 
End If 


xp(n + 1) = xp(n) + (1# / 6#) * (11 + 2# * 12 + 2# * 13 + 14) 
zp(n + 1) = zp(n) + (1# / 6#) * (lzl + 2# * lz2 + 2# * lz3 + lz4) 

xx(n + 1) = xx(n) + (1# / 6#) * (kl + 2# * k2 + 2# * k3 + k4) 
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zz(n + 1) = zz(n) + (1# / 6#) * (kzl + 2# * kz2 + 2# * kz3 + kz4) 

If onoffflag = 1 Then 

D(n + 1) = D(n) + (1# / 6#) * (mdl + 2# * md2 + 2# * md3 + md4) '3/28/03 
Ts(n + 1) = Ts(n) + (1# / 6#) * (mtl + 2# * mt2 + 2# * mt3 + mt4) 

End If 
DoEvents 
If (Halt) Then 
DoEvents 
Halt = False 

Response = MsgBox(Msg, Style, Title) 

If Response = vbNo Then 
evap = 2# 
nl = n 
GoTo 15 
End If 
End If 

If 0 = (n Mod speedup) Or n >= kl Then 
Worksheets("Process").CommandButtonl. Width = ((n / kl) * 165.75) 
Worksheets("Process").CommandButtonl. Caption = ((n / kl) * 100) & "%" 
Worksheets("Process").CommandButtonl. Height = 20.25 
End If 

If flag = 4 And ss = 1 And onoff_flag = 1 Then 
Sheets("Process").Cells(14, 12) = "Estimation done" 

Sheets("Process").Cells(18, 10) = "Temperature is steady at t = " & dtt & "seconds" 

Sheets("Process").Cells(6, 10) = D02 / lambda 

Sheets("Process").Cells(5, 10) = lambda 

dttemp = dt 

dt = dt_old 

'Sheets("Process").Cells(9, 10) = n 
'Sheets("Process"). Cells) 10, 10) = kl 
'Sheets("Process").Cells(7, 10) = dt '9/10/03 
'Sheets("Process").Cells(8, 10) = dt_temp 
GoTo 50 
Else 

If flag = 4 And n >= kl And ss = 0 Then 
Response = MsgBox(Msg_sizel, Style 1, Title 1) 
dt_temp = dt 
dt = dt_old 
GoTo 50 
End If 
End If 

If (D(n + 1) < 0.00000001) And onoff_flag = 1 Then 
evap =1# 
nl = n 
GoTo 15 
End If 
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n = n + 1 

If n <= kl Then 
GoTo 20 
End If 


15 

Sheets("Process").Cells(14, 12) = "Writing data into cells" 

If evap =1# Then 
nf = nl - 1 

Elself evap = 2# Then 
nf = nl - 1 

Elself evap = 0# Then 
nf = kl 
End If 

n = 0 

ttldata = Sheets("Process").Cells(35, 5) 
sd = Sheets("Process").Cells(35, 5) 
speedup = lnt(sd / 10) 

25 

Urm = (((xp(n) - ug) A 2) + ((zp(n) - wg) A 2)) A 0.5 
If onoff_flag = 1 Then Rems = rhog * Urm * D(n) / mu 
If onoff_flag = 0 Then Rems = rhog * Urm * D(0) / mu 

If Rems <= 1000 Then 
If Rems = 0# Then 

Cdm = 0# 

Else 

Cdm = (((Rems A (2# / 3#)) / 6#) + 1#) * 24# / Rems 
End If 
Else 

Cdm = 0.424 
End If 

Sheets("Data").Cells(nnn, 1) = n * dt 
Sheets("Data").Cells(nnn, 2) = xx(n) 

Sheets("Data").Cells(nnn, 3) = zz(n) 

Sheets("Data").Cells(nnn, 4) = xp(n) 

Sheets("Data").Cells(nnn, 5) = zp(n) 

Sheets("Data").Cells(nnn, 6) = Cdm 
If onoff_flag = 1 Then 
Dplot = (D(n) / D(0)) A (2#) 

Sheets("Data").Cells(nnn, 7) = Dplot 
Sheets("Data").Cells(nnn, 8) =Ts(n) 

End If 

If 0 = (nnn Mod speedup) Or n >= nf Then 
Worksheets("Process").CommandButtonl. Width = (((nnn - 1) / ttldata) * 165.75) 
Worksheets("Process").CommandButtonl.Caption = (((nnn - 1) / ttldata) * 100) & "%" 
Worksheets("Process").CommandButtonl. Height = 20.25 
End If 

DoEvents 
If (Halt) Then 
DoEvents 
Halt = False 

Response = MsgBox(Msg, Style, Title) 
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If Response = vbNo Then 
GoTo 10 
End If 
End If 

n = n + userchoice 
nnn = nnn + 1 
If n > nf Then 
GoTo 10 
Else 

GoTo 25 
End If 


10 

last = Second(Time) + Minute(Time) * 60 + Eloiir(Time) * 3600 
Sheets("Process").Cells(16, 10) = last - begin 
If onoff_flag = 1 Then 
Sheets("Evap").Cells(4, 10) = D02 / lambda 
Sheets("Evap").Cells(3, 10) = lambda 
End If 

If (evap =1#) Then 
' nnn = nnn + 1 

' Sheets("Data").Cells(nnn, 1) = (nf + 1#) * dt 
' Sheets("Data").Cells(nnn, 7) = 0# 

Response = MsgBox(Msg3, Style 1, Title) 

End If 


Macro 1 'plot data on chart 

'sd = Sheets("Process").Cells(35, 5) 

'print variables in a file (Tecplot format) 

If Sheets("Process").Cells(46, 9) = "YES" Then 
location = Sheets("Process").Cells(23, 11) 

Sheets("Process").Cells(14, 12) = "Writing data into file : " + location 
Open location For Output As #1 

Print #1, "TITLE = "; Spc(2); "Spray Jet In CrossFlow"; """" 

If onoffflag = 1 Then 

If Sheets("Process").Cells(47, 9) = "millimeter" Then 
unitconv = 1000# 

Print #1, "Variables ="; Spc(2); "Time(sec)"; Spc(2); "X(mm)"; Spc(2); "Z(mm) 

Spc(2); "U(mm/s)"; Spc(2); "W(mm/s)"; Spc(2); "Cd"; Spc(2); 

"D(Normalized)"; Spc(2); "T(K)"; """" 

Else 

unitconv =1# 

Print #1, "Variables ="; Spc(2); "Time(sec)"; Spc(2); "X(meter)"; Spc(2); 

"Z(meter)"; Spc(2); "U(m/s)"; Spc(2); "W(m/s)"; Spc(2); "Cd"; Spc(2); """" 

"D(Normalized)"; Spc(2); "T(K)"; """" 

End If 
Else 

If Sheets("Process").Cells(47, 9) = "millimeter" Then 
unitconv = 1000# 

Print #1, "Variables ="; Spc(2); "Time(sec)"; Spc(2); "X(mm)"; Spc(2); "Z(mm) 

Spc(2); "U(mm/s)"; Spc(2); "W(mm/s)"; Spc(2); "Cd"; """" 

Else 

unitconv =1# 

Print #1, "Variables ="; Spc(2); "Time(sec)"; Spc(2); "X(meter)"; Spc(2); 

"Z(meter)"; """"; Spc(2); "U(m/s)"; Spc(2); "W(m/s)"; Spc(2); "Cd"; """" 
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End If 
End If 
sdd = 0 


'plot data in a file (Tecplot format) 

Print #1, "ZONE I ="; nnn - 2 & Spc(2); "F = POINT" 

nnn = nnn - 1 
For sf= 2 To nnn 
al = Sheets("Data").Cells(sf, 1) 
bl = Sheets("Data").Cells(sf, 2) 
b 1 = b 1 * unitconv 
cl = Sheets("Data").Cells(sf, 3) 
cl = cl * unitconv 
dl = Sheets) "Data" ).Cells(sf, 4) 
d 1 = d 1 * unitconv 
el = Sheets) "Data"). Cells(sf, 5) 
el = el * unitconv 
fl = Sheets)" Data" ) .Cells(sf, 6) 

If onoff_flag = 1 Then 
gl = Sheets("Data").Cells(sf, 7) 
hi = Sheets) "Data"). Cells(sf, 8) 

End If 

If onoff_flag = 1 Then 

Print #1, al; Spc(2); bl; Spc(2); cl; Spc(2); dl; Spc(2); el; Spc(2); fl; Spc(2); gl; Spc(2); hi 
Else 

Print #1, al; Spc(2); bl; Spc(2); cl; Spc(2); dl; Spc(2); el; Spc(2); fl 
End If 
Next sf 
Close #1 

End If 

Sheets("Process"). Cells) 14, 12) = "Completed" 

50 

SuperExit = False 
End Sub 
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Other Subroutines 


Sub driverforrunge(Ts) 

Call transfernumber(Ts) 

Call referencecond(Ts) 

Call air 

Call fuel_thermal_prop(Ts) 

'fuel-air mixture 

kg = yar * ka + yfr * kfv 

cpg = yar * cpa + yfr * cpfv 

BT = cpg * (Tinf - Ts) / lat 

'evaporation contant 

lambda = 8# * Log(l# + BM) * kg / (cpg * rhof) 
End Sub 


Sub transfernumber(Ts) 

If (otherfuel =1) Then 
LoglOT = Log(Ts) / Log(10#) 

LoglOP = Apv + Bpv / Ts + Cpv * LoglOT + Dpv * Ts + Epv * Ts * Ts 
pfs = 10# A (LoglOP) * 0.13333 'convert 'mmHg' to 'Kpa' 
pfs = Application. WorksheetFunction.Min(p, pfs) 

Else 

pfs = Exp(acc - bee / (Ts - 43#)) 'unit='Kpa' 
pfs = Application. WorksheetFunction.Min(p, pfs) 

End If 

yfs = 1# / (1# + (p / pfs - 1#) * (mwg / mwd)) 

BM = yfs / (1# - yfs) 

End Sub 


Sub referencecond(Ts) 

Tref = Ts + (Tinf - Ts) / 3# 
yfr = 2# * yfs / 3# 
yar = 1# - yfr 
End Sub 


Sub air() 

'air: density 

rhog = 355.91 * Tref A (-1.0032) 'kg/m A 3 
'air: viscosity 

mu = 0.00000007 * Tref A 3 - 0.0003 * Tref A 2 + 0.6227 * Tref + 18.761 
mu = mu * 0.0000001 'Ns/m A 2 

'air: prandtl number 
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Pra = 0.0000000000002 * Tref A 4 - 0.000000001 * Tref A 3 + 0.000002 * Tref A 2 - 0.0009 * Tref + 0.8632 


'air: thermal conductivity 

ka = 0.00000000002 * Tref A 3 - 0.00000005 * Tref A 2 + 0.0000969 * Tref + 0.0008289 'W/(mK) 
'air: specific heat 

cpa = -0.00000000005 * Tref A 3 + 0.0000002 * Tref A 2 - 0.00001 * Tref + 1.0041 'kJ/kg 
cpa = cpa * 1000# 'multiply by 1 kilo => J/(kgK) 

End Sub 


Sub fuel thermal prop(Ts) 

If (otherfuel =1) Then 
'fuel: liquid density 
Tdocrit = Ts / Tcrit 

rhof = 1000# * Arho * Brho A ((-1#) * (1# - Tdocrit) A nrho) 'convert g/cm A 3' to 'kg/m A 3 

'fuel: vapor thermal conductivity 

kfv = Akfv + Bkfv * Tref + Ckfv * Tref* Tref 'W/(mK) 

'fuel: vapor specific heat at constant pressure 

cpfv = Acpfv + Bcpfv * Tref + Ccpfv * Tref * Tref + Dcpfv * Tref * Tref * Tref + Ecpfv * Tref * Tref * Tref * 
Tref 

cpfv = cpfv * 1000# / mwd 'convert to J/(kgK) 

'fuel: liquid specific heat at constant pressure 

cpf = Acpf + Bcpf * Ts + Ccpf * Ts * Ts + Dcpf * Ts * Ts * Ts 

cpf = cpf* 1000# / mwd 'J/(kgK) 

'fuel: liquid latent heat of vaporization 
lat = Alat * (1# - Tdocrit) A nlat 
lat = lat * 1000000# / mwd 'J/kg 
Else 

'fuel: liquid density 

rhof = rhofr * (1# - 1.8 * cexp * (Ts - 288.6) - 0.09 * ((Ts - 288.6) / (Tcrit - 288.6)) A 2#) 'kg/m A 3 

'fuel: vapor thermal conductivity 
nkfv = 2# - 0.0372 * (Tref / Tbn) A 2 

kfv = 0.000001 * (13.2 - 0.0313 * (Tbn - 273#)) * (Tref / 273#) A nkfv 'W/(mK) 

'fuel: vapor specific heat at constant pressure 

cpfv = (0.363 + 0.000467 * Tref) * (5# - 0.001 * rhofr) 'kJ/(kgK) 

cpfv = cpfv * 1000# 'J/(kgK) 

'fuel: liquid specific heat at constant pressure 

cpf = (0.76 + 0.00335 * Ts) * (0.001 * rhof) 'kJ/(kgK) 

cpf = cpf* 1000# 'J/(kgK) 

'fuel: liquid latent heat of vaporization 

lat = latbn * ((Tcrit - Ts) / (Tcrit - Tbn)) * (0.38) 'kJ/kg 

lat = lat * 1000# 'J/(kgK) 

End If 
End Sub 
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Sub taperd() 

fl = Sheets("Process").ComboBox2. Value 

If fl = "Option 1" Then flag = 1 

If fl = "Option 2" Then flag = 2 

If fl = "Option 3" Then flag = 3 

If fl = "Evap time estimate" Then flag = 4 

fl = Sheets("Process").ComboBox3. Value 
If fl = "ON" Then onoff_flag = 1 
If fl = "OFF" Then onoff flag = 0 

'Getting the input 

r = Sheets("Process").Cells(5, 3) 
r = r / 1000000# 

rhofr = Sheets("Process").Cells(8, 3) 

Tcrit = Sheets("Process").Cells(9, 3) 

Tbn = Sheetsf "Process"). Cellsf 1 0, 3) 
cexp = Sheets("Process").Cells(ll, 3) 
latbn = Sheets("Process").Cells(12, 3) 
mwd = Sheets("Process").Cells(13, 3) 
acc = Sheets("Process").Cells(14, 3) 
bcc = Sheets("Process").Cells(15, 3) 

xxi = Sheets("Process").Cells(17, 3) 
xpi = Sheets("Process").Cells(18, 3) 
zzi = Sheets("Process").Cells(19, 3) 
zpi = Sheets("Process").Cells(20, 3) 

Tinit = Sheets("Process").Cells(21, 3) 

mwg = Sheets("Process").Cells(23, 3) 
ug = Sheets("Process").Cells(24, 3) 
wg = Sheets("Process").Cells(25, 3) 

Tinf = Sheets("Process").Cells(27, 3) 
p = Sheets("Process").Cells(28, 3) 
g = Sheets("Process").Cells(29, 3) 


dt = Sheets("Process").Cells(31, 3) 
kl = Sheets("Process").Cells(32, 3) 

If flag = 4 And onoff_flag = 1 Then 
dt_old = dt 
dt = 0.000001 
kl = 400000 

If r > 0.00005 Then kl = 800000 
End If 

userchoice = Sheets("Process").Cells(33, 3) 
sstemp = Sheets("Process").Cells(34, 3) 

area noevap = Sheets("Process").Cells(6, 3) 
volume noevap = Sheets("Process").Cells(7, 3) 

End Sub 
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